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ABSTRACT 

We perform 3D lithography simulations by using a finite-element solver. To proof applicability to real 3D 
problems we investigate DUV light propagation through a structure of size 9 /im x 4 /im x 65 nm. On this 
relatively large computational domain we perform rigorous computations (No Hopkins) taking into account a grid 
of 11 x 21 source points with two polarization directions each. We obtain well converged results with an accuracy 
of the diffraction orders of about 1%. The results compare well to experimental aerial imaging results. We 
further investigate the convergence of 3D solutions towards quasi-exact results obtained with different methods. 
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1. INTRODUCTION 

Shrinking feature sizes in optical lithography lead to increasing importance of rigorous simulations for process 
design. 1 Modern lithography simulators include modules describing illumination, transfer of the optical field 
through the mask and aberrating optical system of the lithographic equipment, the propagation inside the 
photoresist, the processes leading to the resist image and - in advanced systems - the etching processes leading 
to the etched image. After nearly two decades of lithography simulation, most of the modules along the simulation 
chain have attained a high degree of maturity. However, the simulation of light propagation through lithography 
masks is still challenging in terms of computational time and memory and accuracy of the results. 

The computation of the print image of a whole chip remains extremely demanding although approximations, 
multi-threading and even hardware accelerators are applied to reduce the runtime of simulations. Rigorous 
simulations are restricted today to small areas and even those simulations suffer from the high computational 
effort. At the same time, the progress on the semiconductor roadmap forces the need of rigorous 3D simulations. 
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Keeping this background in mind, we employed a frequency-domain finite-clement method (FEM) solver 
for Maxwell's equations. In a recent benchmark this solver has been shown to be superior in accuracy and 
computational time requirements by several orders of magnitude, compared to a FDTD solver. 2 Further, 
this solver has been successfully applied to a wide range of 3D electromagnetic field computations including 
left-handed metamaterials in the optical regime, 3,4 photonic crystals, 5 and nearfield-microscopy. 6 



Px 

Figure 1. Schematics of the 3D test structure. The height of the lines is 65.4 nm, the lateral size of the computational 
window is p x x p y — 4 /im x 9 /jm. 

We consider light scattering off a mask which is periodic in the x— and y— directions and is enclosed by 
homogeneous substrate (at z su b) and superstrate (at z sup ) which are infinite in the — , resp. +z— direction. 
However, the presented FEM concept holds as well for non-periodic scattering objects, where the surround- 
ing space is either homogeneous or consists of layered media or waveguide structures. Light propagation in 
the investigated system is governed by Maxwell's equations where vanishing densities of free charges and cur- 
rents are assumed. The dielectric coefficient e(x) and the permeability n(x) of the considered photomasks 
are periodic and complex, e (x) = e (x + a) , \i (x) = fi(x + a). Here a is any elementary vector of the pe- 
riodic lattice. For given primitive lattice vectors d\ and 02 an elementary cell H C K is defined as Q, = 
{x G K 2 I x = OL\d\ + 0:20,2; < ot\, 0L2 < l} x [z S ub, z sup \. A time-harmonic ansatz with frequency u and mag- 
netic field H(x,t) — e~ lut H(x) leads to the following equations for H(x): 

• The wave equation and the divergence condition for the magnetic field: 



V x — - V x H(x) - u 2 fi(x)H(x) 
e(x) 
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• Transparent boundary conditions at the boundaries to the substrate (at z su b) and superstrate (at z sup ), 
dfl, where H m is the incident magnetic field (plane wave in this case), and n is the normal vector on <9f2: 



1 



e(x) 



— V x (H - H m ) x n = DtN(H - H m ), x e <9fl 



(3) 



The DtN operator (Dirichlet-to-Neumann) is realized with the PML method. 7 This is a generalized 
formulation of Sommerfeld's radiation condition; it can be realized alternatively by the Pole condition 
method. 8 
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Figure 2. Aerial image of the test-structure depicted in Fig. 1. Rigorous FEM simulation (JCMharmony) (a), simulation 
using thin mask approximation (Kirchhoff) (b), and experimentally attained image (AIMS) (c) agree well. 



• Periodic boundary conditions for the transverse boundaries, 90, governed by Bloch's theorem 9 : 

H(x) = e ik ' s u(x), u(x) = u(x + a), (4) 
where the Bloch wavevector k G K 3 is defined by the incoming plane wave H m . 



Similar equations are found for the electric field E(x 7 1) 



t E(x); these are treated accordingly. The finite- 



element method solves Eqs. (1) - (4) in their weak form, i.e., in an integral representation. The computational 
domain is discretized with triangular (2D) or tetrahedral/prismatoidal (3D) patches. The use of prismatoidal 
patches is well suited for layered geometries, as in photomask simulations. This also simplifies the geometry 
description of the mask layout. Sidewall angles different from 90deg are not regarded throughout this paper; 
however, they can easily be implemented with reasonable restrictions. The function spaces are discretized 
using Nedelec's edge elements, which are vectorial functions of polynomial order (here, first to fourth order) 
defined on the triangular or tetrahedral patches. 10 In a nutshell, FEM can be explained as expanding the field 
corresponding to the exact solution of Equation (1) in the basis given by these elements. This leads to a large 
sparse matrix equation (algebraic problem). For details on the weak formulation, the choice of Bloch-periodic 
functional spaces, the FEM discretization, and our implementation of the PML method we refer to previous 
works. 5, 7,11 In future implementations performance will further be increased by using higher order elements, 
p > 4, /ip-adaptive methods, and by using elements of different polynomial order parallel and orthogonal to 
the layers of the layered geometry (corresponding to the x — y— plane of Fig. 4 b). To solve the algebraic 
problem on a standard workstation either standard linear algebra decomposition techniques (LU-factorization, 
e.g., package PARDISO 12 ) or iterative and domain decomposition methods 13-15 are used, depending on problem 
size. Domain decomposition methods as shown in Ref. 14 for 2D FEM simulations can be easily transferred to 
layered 3D geometries (typical photomask geometries) and other 3D geometries. Due to the use of multi-grid 
algorithms, the computational time and the memory requirements grow linearly with the number of unknowns. 



2. 3D SIMULATIONS 



We investigate a 3D test structure as schematically depicted in Figure 1. The structure consists of MoSi- lines 
of height h with a sidewall angle of 90deg on a glass substrate. This pattern was chosen to have low vector and 
interference effects but still significant 3D effects, and because it appears in current lithography production. With 



parameter 


data set 1 


data set 2 


Px 


4000 nm 


800 nm 


Py 


9000 nm 


800 nm 


h 


65.4 nm 


65.4nm 


w 


390 nm 


400 nm 


d 


520 nm 




h 


2210 nm 




h 


3910 nm 




h 


6000 nm 




Ao 


193.0 nm 


m 


2.52 + 


0.596i 


n 2 


1.56306 



Table 1. Parameter settings for the simulations in Section 2 (data set 1) and Section 3 (data set 2). 
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Figure 3. Convergence of the zero order far field coefficient for light transition through the test structure described 
in Fig. 1 and Table 1. Magnitude of the far field coefficient in dependence on the number of ansatz functions for the 
numerical solution N for finite elements of polynomial degree p = 1 ... 4. Results of an estimated relative error around 
1-2% are reached for N w 5 ■ 10 5 and p = 4. 

the project parameters as given in Table 1 the size of the computational domain is around 10 3 cubic wavelengths. 
We discretize the computational domain using prismatoidal patches, and we use higher order, vectorial ansatz 
functions as finite elements defined on these patches. 

Figure 3 shows the convergence of the simulated zero order far field coefficient. Plotted is the magnitude 
of the far field coefficient vs. number of degrees of freedom of the finite element expansion. Data points are 
attained using finite elements of first, second, third, and fourth polynomial order and using meshes with different 
refinement levels. With increasing finite clement degree and with decreasing mesh size the results converge. 
From these results we guess the relative error for a solution using fourth order finite elements and using about 
5 x 10 5 unknowns is of the order of 1%. 

To model illumination with a realistic source we construct a grid of 11 x 21 source points in wavevector-space. 
For each source point we define two plane waves with orthogonal polarizations. This makes a total of 462 sources. 
In order to obtain the scattering response of an extended, measured C-Quad source we linearly interpolate the 
results for each measured source point between the closest simulated source points and superpose the results. 



JCMharmony allows to calculate the scattering response of all of the 462 sources in a single programme run. The 
total computation time on a workstation (2 64bit processors, around 20 GB RAM) to obtain the 462 near field 
solutions, each with 482,040 degrees of freedom, was around 100 minutes. We currently work on the significant 
reduction of the near field computation time. 

A simulation on this area is impossible with present Solid-E 3.3 or Prolith 9.1 even for a single illumination 
direction due to memory consumption at necessary resolution. Also the simulation time with unacceptable coarse 
grids is orders of magnitude higher. 

We use the far field coefficients to generate an aerial image. Figure 2 a) shows a pseudo color representation 
of the intensity distribution in the image plane (demagnification factor 4). The scalar intensity distribution has 
been obtained from the vectorial electric field distribution. Figure 2 b) shows a similar intensity distribution 
obtained using a non-rigorous method (Kirchhoff approximation) . Obviously, the approximation compares well 
with the rigorous solution for this specific simulation example. I.e., interference effects, high-NA effects or other 
effects do not play a role here. This example was rather chosen to demonstrate the applicability of a rigoros 
method to large 3D problems. 

Figure 2 c) shows an experimentally obtained aerial image obtained using AIMS. The minima in the experi- 
mental aerial image of the mask structure are more pronounced than in the simulation. This can be, e.g., caused 
by uncertainties in the geometry of the sample. As has been shown in different works 16 the high accuracy and 
speed of rigorous FEM simulations can be utilized to obtain precise informations about the sample geometry or 
material parameters by optimizing the deviation from experimentally obtained data. 

3. VALIDATION OF THE RESULTS 

In order to validate the results of the FEM solver we have performed several tests. Here we present 3D simulations 
of line masks and compare the results to results using 2D simulations of the same physical settings. For these we 
have we have investigated the convergence towards results obtained with various methods, as reported earlier. 2 

Figure 4 shows a schematics of the geometry of a periodic line mask. TM polarized light is incident from the 
substrate and is diffracted into various diffraction orders. The geometry parameters are listed in Table 1 (data 
set 2). Please see Reference 2 for more details on this test case. The 2D simulation results given in the first line 
of Table 2 are the best converged results from this reference for the given test case. We therefore refer to these 
results as 'quasi-exact' results. Table 2 further shows results obtained with 3D FEM on the geometry as depicted 
in Figure 4 b). Vectorial finite elements of order 1 to 4 on grids with different mesh refinements have been used 
to obtain rigorous near field solutions from which the zero order far field coefficients are obtained. As expected 
and as can be seen from the results, most accurate results are obtained using elements of high order and fine 
meshes. Figure 5 shows the convergence of the results on several fixed FEM grids with elements of increasing 
polynomial order. As expected, with increasing polynomial order the numerical approximation error converges 
exponentially towards zero. 

4. CONCLUSIONS 

We have performed rigorous 3D FEM simulations of light transition through a large 3D photomask (size of the 
computational domain about 1000 cubic wavelengths). We have achieved results at high numerical accuracy 
which compare well to experimental findings using aerial imaging. We have checked the convergence for this 3D 
case and we have checked the convergence of the method for a simpler case, where a quasi-exact result is available. 
Our results show that rigorous 3D mask simulations can well be handled at high accuracy and relatively low 
computational cost. 
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Figure 4. Schematics of the computational domain of a periodic linemask for 2D calculations (a) and for 3D calculations 
(b): The geometry consists of a line of width w (at center of the line), height h and sidewall angle f3, on a substrate 
material SiC>2, surrounded by air. The geometry is periodic in ^-direction with a pitch of p x and it is independent on 
the y-coordinate. The refractive indices of the different present materials are denoted by m (line), ri2 (substrate) and U3 
(air), n 3 = 1.0. 
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-0.21943 
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JCMharmony 3D, first order elements (TM) 
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3?(FC ) 


%(FCo) 


\FC \ 
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470 


-0.15665 


O.03066 


0.15963 
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2230 


-O.17802 


O.09517 


0.20186 
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9821 


-0.23282 


0.23871 


0.33345 
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35741 


-0.21775 


0.26485 


0.34287 
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84139 


-O.22011 


0.26071 


0.34120 
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299690 


-0.21847 


0.26942 


0.34687 


JCMharmony 3D, second order elements (TM) 


n 


N 


n(FC ) 


Z(FC ) 


\FC \ 


2 


1884 


-0.26261 


0.22353 


0.34486 


2 


13672 


-0.22516 


0.26312 


0.34631 


2 


54398 


-0.21984 


0.26810 


0.34671 


2 


99468 


-0.21885 


0.27110 


0.34841 


2 


376786 


-O.21910 


0.27145 


0.34884 


JCMharmony 3D, third order elements (TM) 


n 


N 


3?(FC ) 


Z(FC ) 


\FC \ 


3 


5652 


-0.22476 


0. 26408 


0.34678 


3 


27264 


-0.21862 


0.27199 


0.34896 


3 


139659 


-0.21971 


0.27199 


0.34964 


JCMharmony 3D, fourth order elements (TM) 


n 


N 


n(FC ) 


%(FCo) 


\FC \ 


4 


12608 


-0.21726 


0.27308 


0.34896 


4 


100036 


-0.21902 


0.27154 


0.34886 


4 


365384 


-0.21931 


0.27170 


0.34917 



Table 2. Comparison of quasi-exact results (first row, see Reference 2 ) to results obtained using adaptive 3D FEM. The 
3D results converge towards the quasi-exact results for increasing finite element degree n and for increasing number of 
degrees of freedom N (i.e., increasing grid refinement). Real and imaginary parts and magnitudes of the th far field 
coefficients computed with 2D and 3D FEM are given in units [V/m], c/. 2 Converged digits are marked in bold face. 



